function [auxOut,f,errorMes] = numF(params)
 
% Computing the steady 
[auxOut,errorMes] = DSGEmodel_ss(params);
params = auxOut.params;
xss = auxOut.xssTrans;
yss = auxOut.yssTrans;
 
%Unfold params
BETTA= params.BETTA;
B= params.B;
CHI= params.CHI;
CHI0= params.CHI0;
THETA= params.THETA;
DELTA= params.DELTA;
ALFA= params.ALFA;
PHI= params.PHI;
PHIzero= params.PHIzero;
ZI= params.ZI;
KAPAw= params.KAPAw;
ETA= params.ETA;
RHOR= params.RHOR;
PHIpai= params.PHIpai;
PHIpai_1= params.PHIpai_1;
PHIy= params.PHIy;
PHIy_1= params.PHIy_1;
PHIc= params.PHIc;
PHIc_1= params.PHIc_1;
PHIl= params.PHIl;
NU= params.NU;
U0= params.U0;
U0d= params.U0d;
OMEGAz= params.OMEGAz;
OMEGAd= params.OMEGAd;
OMEGAn= params.OMEGAn;
OMEGAp= params.OMEGAp;
OMEGAa= params.OMEGAa;
RHOz= params.RHOz;
RHOd= params.RHOd;
RHOn= params.RHOn;
RHOp= params.RHOp;
RHOa= params.RHOa;
Kss= params.Kss;
OUTPUTss= params.OUTPUTss;
PAIss= params.PAIss;
MUZss= params.MUZss;
Css= params.Css;
AA= params.AA;
lss= params.lss;
Rss= params.Rss;
Wss= params.Wss;
 
% Current values in the model;
% x;
c_ba1 = xss(1);
muz_cu = xss(2);
d_cu = xss(3);
n_cu = xss(4);
paistar_cu = xss(5);
a_cu = xss(6);
% xp;
c_ba1p = xss(1);
muz_cup = xss(2);
d_cup = xss(3);
n_cup = xss(4);
paistar_cup = xss(5);
a_cup = xss(6);
 
% y
c_cu = yss(1);
pai_cu = yss(2);
evf_cu = yss(3);
% yp;
c_cup = yss(1);
pai_cup = yss(2);
evf_cup = yss(3);
 
%% Setting the dimension of f matrices 
f  = zeros(length(xss)+length(yss),1);
 
%% f Function evaluation
% START DISPLAYING f
f(1,1)=(BETTA*Rss*exp(-d_cu)*exp(-pai_cup)*exp(d_cup)*exp(PHIy*log(-(exp(c_cu) + DELTA*Kss)/(OUTPUTss*((ZI*(exp(pai_cu)/PAIss^NU - 1)^2)/2 - 1))) + PHIy_1*log(-(exp(c_cup) + DELTA*Kss)/(OUTPUTss*((ZI*(exp(pai_cup)/PAIss^NU - 1)^2)/2 - 1))) + PHIc*(log(exp(muz_cu)) - log(MUZss) - exp(c_ba1) + exp(c_cu)) + PHIc_1*(log(exp(muz_cup)) - log(MUZss) - exp(c_ba1p) + exp(c_cup)) - PHIpai*(log(exp(paistar_cu)) - log(exp(pai_cu)) + log(PAIss)) - PHIpai_1*(log(exp(paistar_cup)) - log(exp(pai_cup)) + log(PAIss)) + PHIl*log(1/(lss*(-(exp(-a_cu)*(exp(c_cu) + DELTA*Kss))/(Kss^THETA*((ZI*(exp(pai_cu)/PAIss^NU - 1)^2)/2 - 1)))^(1/(THETA - 1)))))*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*(AA/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))*(exp(d_cup)*(((exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))/Css^CHI0)^(1 - CHI)/(CHI - 1) - U0d + (PHIzero*exp(n_cup)*(1 - 1/(-(exp(-a_cup)*(exp(c_cup) + DELTA*Kss))/(Kss^THETA*((ZI*(exp(pai_cup)/PAIss^NU - 1)^2)/2 - 1)))^(1/(THETA - 1)))^(1 - 1/PHI))/(1/PHI - 1)) - U0 + (AA*BETTA)/exp(evf_cup)^(1/(ALFA - 1)))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(2,1)=(Kss^THETA*exp(a_cu)*(THETA - 1)*(ETA + (ZI*exp(pai_cu)*(exp(pai_cu)/PAIss^NU - 1))/PAIss^NU - (BETTA*ZI*exp(-d_cu)*exp(d_cup)*exp(muz_cup)*exp(pai_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cup) + DELTA*Kss)*((ZI*(exp(pai_cu)/PAIss^NU - 1)^2)/2 - 1)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*(exp(pai_cup)/PAIss^NU - 1)*(AA/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))*(exp(d_cup)*(((exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))/Css^CHI0)^(1 - CHI)/(CHI - 1) - U0d + (PHIzero*exp(n_cup)*(1 - 1/(-(exp(-a_cup)*(exp(c_cup) + DELTA*Kss))/(Kss^THETA*((ZI*(exp(pai_cup)/PAIss^NU - 1)^2)/2 - 1)))^(1/(THETA - 1)))^(1 - 1/PHI))/(1/PHI - 1)) - U0 + (AA*BETTA)/exp(evf_cup)^(1/(ALFA - 1)))))^ALFA)/(PAIss^NU*(exp(c_cu) + DELTA*Kss)*((ZI*(exp(pai_cup)/PAIss^NU - 1)^2)/2 - 1)*(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI) - 1))/(ETA*(KAPAw*Wss - (Css^CHI0*PHIzero*exp(n_cu)*((exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))/Css^CHI0)^CHI*(KAPAw - 1))/(1 - 1/(-(exp(-a_cu)*(exp(c_cu) + DELTA*Kss))/(Kss^THETA*((ZI*(exp(pai_cu)/PAIss^NU - 1)^2)/2 - 1)))^(1/(THETA - 1)))^(1/PHI))*(1/(-(exp(-a_cu)*(exp(c_cu) + DELTA*Kss))/(Kss^THETA*((ZI*(exp(pai_cu)/PAIss^NU - 1)^2)/2 - 1)))^(1/(THETA - 1)))^THETA) + 1;
f(3,1)=exp(-evf_cu)*((exp(muz_cup)^((CHI - 1)*(CHI0 - 1))*(exp(d_cup)*(((exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))/Css^CHI0)^(1 - CHI)/(CHI - 1) - U0d + (PHIzero*exp(n_cup)*(1 - 1/(-(exp(-a_cup)*(exp(c_cup) + DELTA*Kss))/(Kss^THETA*((ZI*(exp(pai_cup)/PAIss^NU - 1)^2)/2 - 1)))^(1/(THETA - 1)))^(1 - 1/PHI))/(1/PHI - 1)) - U0 + (AA*BETTA)/exp(evf_cup)^(1/(ALFA - 1))))/AA)^(1 - ALFA) - 1;
f(4,1)=exp(-c_ba1p)*exp(c_cu) - 1;
f(5,1)=RHOz*log(exp(muz_cu)/MUZss) - log(exp(muz_cup)/MUZss);
f(6,1)=RHOd*log(exp(d_cu)) - log(exp(d_cup));
f(7,1)=RHOn*log(exp(n_cu)) - log(exp(n_cup));
f(8,1)=RHOp*log(exp(paistar_cu)) - log(exp(paistar_cup));
f(9,1)=RHOa*log(exp(a_cu)) - log(exp(a_cup));
% END DISPLAYING f
 
end
